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We study the synchronisation properties of the Kuramoto model of coupled phase oscillators on 
a general network. Here we distinguish the ability of such a system to self-synchronise from the 
stability of this behaviour. While self-synchronisation is a consequence of genuine non-perturbative 
dynamics, the stability in dynamical systems is usually accessible by fluctuations about a fixed 
point, here taken to be the synchronised solution. We examine this problem in terms of modes of the 
graph Laplacian, by which the absolute Lyapunov stability of the synchronised fixed point is readily 
demonstrated. Departures from stability are seen to arise at the next order in fluctuations where the 
dynamical equations resemble those for species population models, the logistic and Lotka-Volterra 
equations. Methods from these systems are exploited to analytically derive new critical couplings 
signalling deviation from classical stability. We observe in some cases an intermediate regime of 
behaviour, between incoherence and synchronisation, where system wide periodic behaviours are 
exhibited. We discuss these results in light of simulations. 

PACS numbers: 89.75.Fb,05.45.Xt,89.75.Hc 



I. INTRODUCTION 

The study of the ability of systems to achieve self- 
synchronisation - coherent behaviour purely through the 
mutual interactions of many components - is a major 
theme in complex systems science. The phenomenon 
arises in many contexts, from ecology, biology, physics 
to social systems. In our case the understanding of self- 
synchronisation through complex networks in command 
and control organisations is a major goal [1]. Notable 
early work on collective synchronisation includes that of 
Wiener pj and Winfree [3^. Kuramoto's elegantly simple 
model upon which we focus, encapsulates much of the 
richness expected in self-synchronising systems. Strogatz 
[5] gives a short elegant introduction while Acebron et al. 
[B] provide a more thorough review. 

The Kuramoto model involves one-dimensional phase 
oscillators described by angles 9i coupled on a network 
described by an undirected graph G of size iV and is given 
by the first order differential equation 

= tj, - — ^ Aij sm{e, -Oj). (1) 
i 

Here Aij is the adjacency matrix encoding the graph 
structure, uji are intrinsic frequencies of the oscillators, 
usually drawn from some statistical distribution and K 
is a uniform coupling controlling the degree to which ad- 
jacent oscillators must mutually adjust with respect to 
each other. At zero coupling the system obviously be- 
haves incoherently with each oscillator moving in simple 
harmonic motion according to its intrinsic frequency. In- 
creasing the coupling enables the population to eventu- 
ally break into synchrony, with some of the oscillators 
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locking into a single core that moves with harmonic mo- 
tion according to the frequency average over the whole 
oscillator ensemble. For the case of an infinite, complete 
graph and symmetric and unimodal distribution of in- 
trinsic frequencies Kuramoto was able to derive an exact 
analytic expression for the critical coupling Kf. at which 
synchrony breaks out 



where g(0) is the central value of the distribution g from 
which the intrinsic frequencies of the individual oscilla- 
tors are selected. 

On a general network no equivalent analytic expression 
for Kuramoto's critical coupling has been derived. This 
paper shall take some steps in that direction. While an- 
alytic methods have been applied to mean-field versions 
of the model, such as in [7], numerical simulations have 
delivered insights into the pattern of synchronisation for 
the complete form of the Kuramoto model on finite gen- 
eral networks. For random graphs (such as Erdos-Renyi) 
it is observed by ^ that, as the coupling is increased, 
synchrony first develops locally on separate hubs which 
finally merge into a single locked system after the cou- 
pling has crossed a final threshold. The same authors find 
that for "scale-free" networks [S] the synchrony develops 
at one point in the graph and then extends to the whole 
graph as coupling is increased. For the Watts-Strogatz 
"small world" network [10] synchronisation occurs even 
for small rewiring probability Such work already 

shows that quite specific graph substructures play a sub- 
tle role in the transient towards synchronisation, both 
in manifesting local synchronisation and, possibly, stim- 
ulating it in other parts of the network. Further insight 
is offered by Arenas et al. [12] who find evidence that 
the relevant sub-structures are the modes of the graph 
Laplacian, an object that is ubiquitous in network cou- 
pled dynamical systems [T3J[TJ. Indeed, for graphs with 
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hierarchical community structure the order in which sub- 
structures settle into synchronisation is according to the 
Laplacian eigenvalue, with the largest mode synchronis- 
ing first and the smallest stabilising last \T^. A more 
thorough review of numerical studies is provided by |15j . 

At the heart of it, synchronisation and its stability 
are two separate phenomena in such systems. These 
aspects are often confused in the literature, so that, in 
some cases, stability criteria (derived perturbatively) are 
treated as synonymous with synchronisability (a deeply 
non-perturbative phenomenon). Intuitively, the distinc- 
tion can be appreciated from the perspective of a poten- 
tial energy landscape for the entire system. The evolution 
of the system from some random initial condition to com- 
plete synchronisation can be visualised as a walk through 
this landscape, from the "edges" (incoherent behaviour) 
to the "centre" (coherence) through a typically rough in- 
termediate landscape. Synchronisation is the means by 
which the system passes through the bulk of the rough 
landscape close to the centre. Stability is a local prop- 
erty of the vicinity of the centre. In the vicinity of the 
synchronised solution, for stability, perturbative methods 
are applicable. This is where this paper will focus and 
provide new analytical insights. The mechanism for the 
system finding its way through the landscape is a genuine 
nonperturbative problem for which analytic methods are 
lacking in the finite general network case though numeri- 
cal studies have provided significant insights. This paper 
will offer an hypothesis, building on the insights from the 
existing numerical studies. 

While most stability studies are performed for the in- 
finite network case using the formalism of probability 
densities [T6|, in the finite case there is little that has 
been published specifically for the Kuramoto model. We 
suspect that this is because many researchers may have 
implicitly appreciated that the synchronised solution is 
Lyapunov stable and decided that there was little more 
to be said. We shall review this aspect in terms of eigen- 
modes of the graph Laplacian. In these terms, the syn- 
chronised solution turns out to be just the zero mode 
of the Laplacian while the normal mode fiuctuations, to 
first order, exponentially decay to zero consistent with 
the Lyapunov stability criterion. More generally, in this 
framework, we make the observation that the synchro- 
nisation problem can be reframed as a population com- 
petition dynamics. The question then becomes: how do 
the normal mode interactions lead to an "extinction" of 
these modes? We answer this question by extending con- 
siderations to second order in the lowest lying Laplacian 
normal mode fiuctuations and extracting equations more 
commonly seen in species-population models. Though 
we argue that all the low-lying Laplacian normal modes 
are significant in the generation of instabilities in the sys- 
tem, we shall study analytically the role, respectively, of 
one and two of the lowest-lying Laplacian normal modes 
in generating instabilities of the synchronised fixed point. 
The lowest Laplacian normal mode, whose relationship to 
measuring connectivity in a network was noted by Fiedler 



[T7], has often been argued to be relevant to the problem, 
as we shall review below. In our analysis we show that 
the dynamics of a single lowest mode are governed by a 
modification of the logistic equation while for two lowest 
modes a form of Lotka- Volterra equations are obtained. 
Though, for many complex networks, these are rather ar- 
tificial cases, we derive critical couplings here to illustrate 
the types of dynamics that can occur. Indeed, we identify 
in some cases a previously (to our knowledge) unnoticed 
regime lying between incoherence and coherence, a type 
of "edge-of-chaos" patterned dynamics. We provide sup- 
port for this by showing examples of simulations for this 
system. 

The paper is structured as follows. In the next section 
we recast the Kuramoto model in terms of modes of the 
graph Laplacian and review the Lyapunov stability of the 
system. Then we develop the problem to second order in 
fiuctuations and study several cases for the structure of 
the lowest normal modes of the Laplacian. We derive 
a number of critical couplings for several special cases. 
We subsequently analyse the behaviour of an order pa- 
rameter for synchronisation of the system in these terms 
and discuss this in light of simulations. The final section 
will discuss our insights into the general mechanism for 
synchronisation and the subtleties of generalising these 
results to the complete infinite graph, originally studied 
by Kuramoto. 

II. LAPLACIAN EIGENMODES AND 
STABILITY 

A. Comparison of Kuramoto model with other 
coupled dynamical models 

Let us begin by contrasting the Kuramoto model, 
Eq.([T|), with another class of coupled dynamical system 
which has been studied intensively for some time: 

^{<|^^) + cJY,A,,[g{^,)^g{(^,)]. 

J 

We shall consistently use the graph (or combinatorial) 
Laplacian [T3j[T4], 

Lij — diSij -^ij 1 (^) 

with di is the degree of node i. The spectrum of the 
Laplacian is positive semi-definite, so that eigenvalues 
can be ordered from smallest to largest: 

= Ao < Ai < ••• < Aat^i. (4) 

(In some of the literature the numbering of modes starts 
at one and not zero.) In these terms, the dynamical 
system can be rewritten 

(i^i^ f{(l)z) - cr^L^jg{(j)j). (5) 
j 
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Synchronisation is represented by fixed point solutions 
4>i = s oi Eq. ([5| . Fluctuations rji about s satisfy then 



Going to a basis of eigenmodes for the fluctuations leads 
to a criterion for stability of the synchronised solution: 



dT[f'{s{T))-ag'{s{T))K]<Q 



for all Laplacian modes r. This "master stability equa- 
tion" was derived by Pecora and Carroll [18 . Many oth- 
ers have exploited this approach to study a variety of 
networks, for example [19 . Significantly, we see that, de- 
pending on the network and the nature of the functions / 
and (7, there may be directions in Laplacian mode space 
for which this condition is not satisfied. Therefore, within 
a first order perturbation calculation we can determine a 
condition for which synchronisation can exist and is sta- 
ble. It is tempting on this basis to see synchronisation 
purely through the lens of Lyapunov stability. 

Unfortunately, the Kuramoto model Eq. ([ij is very dif- 
ferent from this system. This is because the interaction 
in the Kuramoto model is a Junction of differences rather 
than a difference of functions. Practically, it means that 
the dynamics are not automatically diagonalised in the 
Laplacian eigenmode basis. However, small fluctuations 
about the synchronisation point can be diagonalised in 
terms of the Laplacian. 



B. Representation in Laplacian Basis 

Before diving into this, let us present some more de- 
tails on the Laplacian. The degeneracy of zero eigen- 
modes of the Laplacian corresponds to the number of 
disconnected components of the graph. We shall concern 
ourselves with connected graphs, so there will be a single 
zero mode. The corresponding zero mode eigenvector, 
within a system of orthonormal eigenvectors, is 



iO)T 



1 



(1,1,. ..,1), 



(6) 



The first nonzero eigenmode, here corresponding to Ai, 
is called the Fiedler with the eigenvalue itself known as 
the algebraic connectivity [17] . The eigenmode itself is 
used in other contexts as a measure of network fragility 
[20| . Though we are dealing with unoriented graphs, it is 
convenient to rewrite the Laplacian in terms of a signed 
incidence matrix B as follows. Let links, which may be 
assigned an arbitrary direction for this purpose, be la- 
belled by the index a. Matrix elements of B are Bia 
where 



according to whether link a is incoming to, not connected 
with or outgoing from node i. The Laplacian is then 
given by 

L = BB'^. 

The signed incidence matrix has deep mathematical 
properties, but one role that it plays is to transform be- 
tween objects respectively in the node and link spaces, 
including the eigenvectors: 



(7) 



In terms of the original eigenvectors we can expand a 
phase angle of the Kuramoto model in the Laplacian basis 



r=0 



Inserting this in the equation of motion Eq.Q we arrive 
at equations for the components: 



ao = VNCot (8) 
dr{t) = c.W-a5]e('-)sin(^e(^)a,(t)); 

a \s>0 ) 

r ^ (9) 



where from this point we summarise the coupling as ct = 
K/N . Also implicit in Eq.([9| is that we treat the uji as 
an A'^-dimensional vector, with Co the average over the 
distribution and 



Bi, 



+1,0,-1 



(r) (r) 



represents the projection of the vector of frequencies Ui 
onto the r-th Laplacian mode. 



C. Interpretation of normal mode dynamics 

From Eq.(|8|, we see that the zero mode is just the syn- 
chronised mode and that it completely decouples from the 
normal modes. Therefore, in terms of Laplacian modes 
there is always a component moving harmonically with 
the average frequency. Posed this way, the concern ceases 
to be how synchronisation is generated but how the nor- 
mal mode dynamics in Eq.Q are suppressed. Indeed, we 
can reformulate the problem in the language of another 
classical complex systems problem: that of inter/intra- 
species population dynamics. Recall that in this context, 
populations can grow, die out and compete with others 
in predator-prey relationships. Within such interactions, 
oscillating population numbers for co-existing species can 
occur. 

In this language, then, pure oscillator synchronisation 
is identical to normal eigenmode population extinction 
and complete oscillator incoherence is identical to mul- 
tiple species normal eigenmode population co-existence. 
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Later in the paper, we shall refine this relationship fur- 
ther and identify specific cases of behaviours according to 
interacting population dynamics manifested here in the 
competing Laplacian normal eigcnmodes. 



D. Lyapunov Stability 



Jadbabaie et al. (5T] using different methods. These au- 
thors also derive certain bounds on the coupling, related 
to our bound Eq.( 13 ), which, however, do not readily re- 



turn the special case derived by Kuramoto for the infinite 
complete graph. We suspect this is a consequence of the 
absolute Lyapunov stability of the system and return to 
this at the end of the next section. 



We can now review standard stability arguments. Let 
the system be close to complete synchronisation, now 
characterised by the zero mode, 

Oiit) « ao{t)vl°'' + small, 

with small fluctuations being some superposition of the 
normal modes. Expanding the interaction term for small 
normal modes to leading order, we obtain 



(r) 



aXrarit); r ^ 0. 



The solution to this is just 



ar{t) 



(1 



)+a,(0)e-" 



(10) 



(11) 



Part of this equation was obtained by Arenas et al. |12| , 
who only presented the second term as their focus was on 
the hierarchy of suppression of modes. This supports the 
idea that the sequence of clusters achieving synchroni- 
sation occurs sequentially in time through the hierarchy 
of eigenvalues in Eq.Q. Specifically, fiuctuations in the 
highest mode of the Laplacian, corresponding to A^v-i, 
are suppressed most quickly; the last mode to be sup- 
pressed is the lowest non-zero mode, the Fiedler with 
eigenvalue Ai. //suppression is achieved, the modes be- 
come constant 



v(°°) - 



fj Ar 



(12) 



Note that the complete graph has a completely degen- 
erate Laplacian normal mode spectrum: \i = N, i = 
1, . . . , N — 1. The implications of this degeneracy are dis- 
cussed in the final section. We observe in Eq.(12) that 



the regime of applicability of synchronisation is charac 
terised by small ai°°^ , namely 



a > 



At- 



(13) 



III. SECOND ORDER APPROXIMATIONS 

In the following we extend the above approximation 
to next order in the Taylor expansion of the interaction 
in the equations of motion and making certain simplify- 
ing assumptions on the structure of the lowest modes of 
the Laplacian. Specifically, we expand the sine in Eq.([9]) 
out to the cubic term, sinx = x — This means 

that our approximations are valid for values of ar larger 
than considered in standard stability analysis thus far, 
out to values ofa;'--^7r/2>l. Furthermore we shall as- 
sume sufficient time to have elapsed in the running of the 
system dynamics so that most Laplacian modes are sup- 
pressed to ai°°\ The exception to this will be a number 
of the lowest modes which will be taken to be dynamical. 
In other words, we exploit the observation in numerical 
simulations that modes with small Laplacian eigenvalue 
synchronise later, though we make no assumptions about 
the order of synchronisation of the higher modes. We 
shall now explore the impact of these few small eigen- 
value dynamical modes on stability. 



A. One dynamical mode 

1. Exact solutions 

In this case we have a single mode that remains time- 
dependent, namely that corresponding to the Fiedler, 
ai{t). All other modes r > 1 assume their steady-state 
values Eq.(12). We assume in this part that there is no 



degeneracy for Ai, which will be relaxed in the next sec- 
tion. Thus we separate the mode sum in Eq.([9]) into two 
parts. 



which can be seen as a (very) weak bound on critical 
values of the coupling such that synchronisation can oc- 
cur. However, this bound does not indicate what types 
of transitions can occur. 

We see that in the vicinity of the synchronized solution 
the system is Lyapunov stable in all directions in Lapla- 
cian mode space: the Lyapunov exponents are exactly the 
Laplacian eigenvalues multiplied by the coupling so that 
there are no regimes of coupling or structure whereby 
they may change sign. This has also been observed by 



s>0 



=(1) 



ai{t) 



S>1 



(14) 



and insert this into Eq. (|9| . Within the expansion of the 
sine we drop all terms cubic in any given mode and those 



quadratic in the constant modes a, 
we extract the equation 



(oo) 



To this order then 



di{t) = w'"'"'' — crAiQ;i(i) -I- axiai^t)"^ 



(15) 
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where we define constants summarising the dynamical 
and structural content of the system 



S>1 

E 



This reflects a genuine point at which the dynamics 
changes distinctly in character in contrast to the bound 
Eq.( 13 ) derived purely within a first order stability anal- 
ysis. 



2. Equilibrium analysis 



Eq.( 15 1 is a variation on the well-known logistic equation. 



We further introduce the notation (the first being the 
discriminant of a quadratic equation) : 



A(0) 



2(7Xiai(0) — (crAi -t- 



/AT) 



2axiai{0) — (crAi 



/AT) 



Even though we have analytical solutions here it is 
worth reinforcing these results by looking at the equilib- 
rium dynamics in this regime. Setting di = in Eq.( 15 1, 
two solutions result: 

± ctAi t vat 

«i — ^ • 

In fact these match on to the analytical solution for t 



The solution to Eq.(15l depends on the sign of the dis- ±00 of Eq.(16) 



criminant Ai. When the coupling strength a is varied, 
say, from large to down small values the discriminant 
changes sign, from positive, through zero, zero to nega- 
tive values. 

For large cAi the solution is 



Ml + VAT) - (aAi - v^)Pi(0)e- 
2(7X1 (1 - Pi(0)e^*) 



(16) 



This will return the first order solution when Ai <C A2. 
Then xi ~ and the right hand side of Eq.(16l can 
be expanded for small xi giving Eq. ( [TT|) . This refiects 
the behaviour that the solution to Eq.(|15[) starts from 
the initial condition and then is exponentially suppressed 
to the constant value Eq.(12). There can occur points 
where, at some fixed time, the denominator of Eq.(16) 
vanishes; the solution shoots to negative infinity and then 
returns to stabilise at the constant value Eq.(12). 



When the discriminant is negative the solution takes 
the form 



ai(<) 



1 



2axi 



aXi + tan ^arctan(2(Ta;iQ;i(0) — crAi) 



(17) 



This displays both "accidental" divergent behaviour, 
when the argument of the tangent becomes tt/2 and, 
more significantly, periodicity in time. In this regime 
then, the normal modes are dominated by a limit cycle 
and the lowest Laplacian mode never converges to steady 
state. 

Therefore the condition Ai = separates two distinct 
regimes of behaviour: for Ai < complete synchroni- 
sation is not attained while for Ai > convergence to 
the fixed point occurs. We can thus extract a critical 
coupling: 



S>1 



AiAg 



(18) 



= lim ai{t). 



As a preview to the next section, let us examine stability 
about these equilibria. Inserting 



"1 



X 



into Eq.(15l we obtain the order x equation 
X = {2axia^ - crAi)x'^ 

= tVaTx* 

The solution manifests a Lyapunov-like exponent: 

X± = X*(0)e^^* 

Again, we are only interested in the t 00 behaviour, 
namely in x^- Examining the strong coupling case a > 
1 /Ai we obtain 



ctAi 

Ai _ ^ 
xi aXi 



(19) 
(20) 



Thus, for t — > cx), Eq.(19l coincides with the original 



lowest order steady state result, Eq.(12) 



We now appreciate the two regimes of behaviour seen 
in Eqs.( 16117 ) in another light: Ai > manifests sta- 
bility, namely convergence as t —^ 00 to the equilib- 
rium fixed point while Ai < gives an imaginary phase, 
namely a limit cycle. 



B. Two dynamical modes 

Let us consider initially now the lowest two normal 
modes, ai and 02, to be dynamical with all others having 
reached steady state. Initially we consider all terms linear 
in the constant modes al°°'^ but up to quadratic in the 
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dynamical modes as{t) (s = 1,2). The full system to this 
order has the form 

di(t) = oj^^'^ - aXiai{t) + (TalXi+ aai{t)a2{t)Yi 



a2 



0-02X2 + aai{t)a2{t)Y2 



-aai it)'^ Z2 



(21) 



with a set of dynamical and structural constants 



v\ J- 


s>2 " 


Yr 


s>2 * 




s>2 * 




a 




a 


Els 


- E-i'nei^^)^ei^) 


D2s 


a 

= Els 


E2s 





by analogy with Xr and c^s for the one-mode case. 
Eqs.(21) are not amenable to analytic solution. How- 



ever, analogously to the analysis of Sect|TTTX2] the equi- 
libria and their stability are open to analysis (as is stan- 
dard in predator-prey models). In other words, the in- 
tention is to stud y th e case where = 0. Denoting such 
solutions to Eqs.(21l by ai"\ where u represents some 



label for the different possible solutions, we then exam- 
ine fiuctations Xr"' about these, namely solutions of the 
form Ur = ai"^ +Xr"'^ keeping only terms of order x- The 



system Eqs.(21 1 with LHS set to zero is soluble however 
the solutions are extremely lengthy and cataloging the 
type of behaviours that can occur is difficult, though we 
should expect the full range of dynamical possibilities: 
stability, limit cycles, tori and instability. To illustrate 
this we consider a further truncation of the system, where 
only cross-interactions are kept. We comment later on 
the possible generality of our insights from this analysis. 



1. Cross-interactions only, non- degenerate case 

Treating the case that Ai ^ A2, in the first instance, 
we are down to the coupled equations: 

di{t) = uj'^^'> - a\iai{t) + aai{t)a2{t)Yi 

d2it) = uj'-^^ ~ a\2ai{t) + aai{t)a2{t)Y2, (22) 



Apart from the first term, we recognise in Eqs.(22) 
the two dimensional Lotka-Volterra system for predator- 
prey dynamics [22^ . That we have arrived at the next 



step in generalisation of population models after the lo- 
gistic equation is consistent with our picture of the nor- 
mal mode dynamics. 



The solutions of Eqs.(22l with di = d2 = are com- 
pactly represented by introducing a new discriminant: 

A2 = (o-AiAa - {uj^^^Y2 - w^^^Yi))^ - Aa\iX2UJ^^^Yi. 



Eqs.(22l exhibit two equilibria: 

(TA1A2 + (^(1)^2 - LJ^^^Yi) =F ■ 



„(±) 



„(±) 



crAiA2 



2ctAiX2 



2CTA2Y1 



(23) 



(u) 

We consider now the fluctuations Xr about these so- 



lutions by insertion into Eq.(22). We arrive at a fluctua- 



tion matrix whose eigenvalues give Lyapunov exponents 
for the equilibria Eqs.(23) which have the form 



'1 



(±) _ 



'2 



with 



4A 

1 

4A1A2 



A = (Ai + A2)VA^ 



B = 



^(A±WtCt{A±B)) (24) 



o-AiA2(Ai + A2) 

-(A2-Ai)(a;«r2-^<')n) 



(25) 



C 



To get a handle on the equilibria in the absence of 
explicit analytic solutions it is worth examining flrst the 
strong coupling limit cr > 1/Xr- For the equilibria we 
obtain from Eqs.(23l 



,.(+) 



„(+) 



„(-) 



„(-) 







aXi 








crA2 




A2 




Y2 


<jXiY2 


Ai 




Yi 


'7X2Y1 



while the Lyapunov exponents Eqs.(24l give 



(+) 



;(+) 



4- 



— (tAi 
-(TA2 
—a \J A1A2 



(-) 



+a\JXxX2 



We recognise then that 01^^ converge to the t ^ oo form 
of the lowest order solutions. Our concern being with 
t > we therefore discard the solutions olt 
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Noting that A,C > 0, we test the possible signs of the 
exponents in Eqs.(24l. Consider, firstly, the case of the 
discriminant A2 > 0. Then it is straightforward to derive 
the following cases: 

A + B>0 < 0> 4^^ < 

A + B < > or complex, l'^'^^ > or complex. 

The complex values occur when {A + B)'^ < C. There- 
fore, for A2 > 0, Qfi and a2 will either be both stable or 
both display limit cycles. In the latter case, there can 



be real parts with positive sign in 4^'' of Eqs.(|24|, indi 



eating instability. We shall illustrate this in the simpler 
degenerate case below. 

More interestingly, there is a precise coincidence be- 
tween the "critical point", given by A2 — 0, and an 
"edge-of-chaos" regime defined by vanishing Lyapunov 
exponents, ^2^'' = 4^"* = 0. The boundaries of this regime 
are governed by the two solutions to A 2 = 0: 

1 



A1A2 

Evidently (t_ < cr+ and the two values indicate sign 
changes in A2: positive for a < a- and a > cr+ and 
negative for (t_ < cr < ct+ . Noting that at cr = we have 
A2 = (w(^)y2 - w(2)y^)2 > then either both critical 
points, a±, are positive or both negative (and therefore 
unphysical as we consider only locally attractive interac- 
tions) : if one critical point appears then there will be a 
second. 

We conclude therefore that between incoherence for 
(7 « and coherence for cr ^ l/A^ there will be an in- 
termediate regime characterised by oscillatory behaviour 
whose frequency is not related to any individual oscilla- 
tor but to the collective behaviour of the low-lying graph 
Laplacian modes. 



2. Degenerate case 

In the case that two lowest Laplacian eigenvalues are 
equal, Ai = A2, a further simplification occurs. The Lya- 
punov exponents relevant for t > turn out to be: 



(+) 
1 

(+) 




A2-crA^|-|-(VA2-|-crA?) 



Consider firstly the case of A 2 > 0. For \/A2 > crA^ 
we have: 



When the discriminant is negative we must take care 
with the modulus. Writing 



^2 = i 



for this case, we have 



/A2-<tA?| = \<jXI - i^/\A^\ 



A2| + a2At. 



Thus 



'1 — 



1^+^ - 



1 
2X1 
1 

2X1^ 



A2\+a^\t + a\l + i^\A2\ 



A2\+a^Xt-c7Xl-i^\A2\ 



The fixed points are therefore "hype rbolic" , havi ng non- 
vanishing real parts. Indeed, since •\/|A2| -I- a'^Xf — aXf > 
for non- vanishing discriminant , the sign of the real part 
of the second mode will evidently be positive signalling 
a genuine instability. 

In summary then, for A2 > we obtain stability, 
namely convergence to the steady state while for A2 < 0, 
we obtain limit cycle behaviour for both modes mixed 
with instability in the second mode. Once again, A2 = 
separates two distinct regimes of behaviour, giving two 
possible critical couplings: 



^{J^^Y2 + Lu'^^'^Yi ± 2y/u;Wu;(^)YiY2). 



A? 



The above analysis can be reproduced here. 

Finding manageable explicit solutions beyond these 
truncations is difficult. However, at a general level we 
can identify the structures which are significant in the 
results we have seen. Typically the critical couplings 
arise, in the cases considered here, from the discriminants 
corresponding to quadratic equations for the equilibria. 
The equilibria for the general system Eqs.(21 1 evidently 



will satisfy higher degree polynomial equations, for which 
corresponding discriminants occur. For example, treat- 
ing the first of Eqs.(21 1 as an equation for ai in terms of 
a2 we have initially a quadratic equation. Substitution 
of ai into the second equation will lead to an equation 
involving square roots of the discriminant which can be 
eliminated by squaring the equation. This will gener- 
ate a quartic with its own discriminant. We thus fore- 
see cases where the Lypunov exponents for the various 
equilibria will develop imaginary parts when the discrimi- 
nants are less than zero. What is more difficult to oversee 
is whether there are cases when the Lypunov's are purely 
real and positive, signals of actual chaotic instability to 
this order. 



(+) 



1 

(+) 



Ai 

4+^ -(tAi. 
On the other hand, for \/ A2 < crX^ the exponents swap. 



C. More dynamical modes 

Going beyond the case of two dynamical modes at 
present appears difficult within analytical approaches. 
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However, the identification of the dynamics as being 
essentially a modification of multi-species competitive 
Lotka-Volterra equations enables us to exploit a general 
result by Smale . Here it is known that for systems of 
more than five population species the equilibria can be 
of any type: fixed-point, limit cycle, torus or attractors. 
However, this is just a statement that anything is possi- 
ble. Thus, whereas thus far we have only seen limit cy- 
cles prohibit complete synchronisation, for more than two 
closely spaced low-lying Laplacian modes we can expect 
tori and genuine chaos to be exhibited in the dynamics. 
The more severe these behaviours, the further from par- 
tial synchronisation and the close to true incoherence the 
system approaches. We can also expect subtle dynamical 
effects, as seen in multi-species Lotka-Volterra systems, 
such as purely interaction generated spatio-temporal pat- 
terns (or spontaneous symmetry breaking) |24j . What we 
lack at this stage are criteria that may give insights into 
the boundaries between different regimes. 

Also, our considerations cannot yet be extended to the 
case of the complete graph: the above approximations 
assume a separation between static and dynamic modes 
which cannot be reconciled with the total degeneracy 
of the Laplacian eigenvalue spectrum for the complete 
graph. The insight that the normal mode dynamics can 
be mapped onto population models is still valid neverthe- 
less, though the specific modified Lotka-Volterra equa- 
tions considered above cannot be taken to hold. For that 
reason it is premature to seek to extrapolate the above 
results for critical couplings to check against Kuramoto's 
analytic result. 



IV. ORDER PARAMETER 

We now explore the impact of the above solutions on 
an order parameter characterising the possible phase of 
the system. Jadbabaie et al. [21] propose a generalisation 
of Kuramoto's order parameter for the general network: 



order parameter is constant 



iV2 



with E the number of edges in the network. This can 
be rendered more transparently by using the identity 
cosx = 1 — 2sin^a;/2, converting to the Laplacian basis 
and explicitly using the form of the zero mode eigenvector 
(which makes it drop out of r^), giving 



(26) 



We thus see explicitly that complete synchronisation as 
characterised by complete extinction of normal modes, 



Qfr = 0, gives r 

tion than the steady state behaviour, ar^' , for which the 



1. However this is a stronger condi- 

(oo) 



(r 



(oo) 



4 

iV2 



However, note that in all our approximations thus far we 
have only taken up to (ctr)'^ to be non-negligible. We 
therefore need only consider, at most, the leading term 
in the sine, due to it being squared in r^. Hence, 



a \ r^O j r 



r'^Q 



where only modes r which are dynamical contribute to 
this. For the steady state case where all modes are static 
no normal modes contribute: 



(■y,(00)-)2 



1 



and steady state is consistent with complete synchroni- 
sation. 

For a single dynamical mode we have 



1 



1 



while for two modes the result is: 



1 



{Xial{t) + \2al{t)). 



(27) 



(28) 



The generalisation to many dynamical modes is self- 
evident, though our analytic methods do not help at this 
stage. 

Drawing on our solutions to the various cases above 
we obtain either « 1 or signals of periodic behaviours 
in the order parameter. Therefore, if the system reaches 
a state where all but the two lowest Laplacian normal 
modes are suppressed a rich range of behaviours can oc- 
cur in the future evolution of the system: complete syn- 
chronisation or partial synchronisation with the modes 
settling into limit cycles. Critical values of the coupling 
can be determined signalling the boundary between os- 
cillatory behaviour and complete synchronisation. The 
approximation underlying these calculations is consistent 
with a ~ 1 — 2 so that these effects should be detectable 
for either large graphs with very high algebraic connectiv- 
ity, Al >> 1, or small graphs with algebraic connectivity 
of order unity. 



V. SIMULATION STUDIES 

The aim of this section is not to show detailed sim- 
ulation studies testing the particular form of the criti- 
cal coupling; our approximations do not warrant at this 
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stage such a study. Rather, the intent is to demonstrate 
cases of the transition between incoherence, oscillatory 
behaviour and synchronisation consistent with the qual- 
itative picture we have developed. 

We use NetLogo 3.1.3 [25], a multi-agent modelling en- 
vironment developed by Uri Wilensky (a version 4.0.4 is 
currently available). In particular, this tool provides a 
number of models for the study of complex networks as 
well as for synchronisation of fireflies. Specifically, we 
use an extension of one of these models developed by 
Dekker which exploits the Kawachi [27] process for 
generating a spectrum of well-studied complex networks 
by variation of a single parameter p, a rewiring probabil- 
ity, starting with a regular graph. For p ~ the graph 
is a ring, for p = 0.05 — 0.1 small world networks are 
produced, for p = 1 — 2 the graph is random (not quite 
Erdos-Renyi) and for p ~ 5 the graph is scale-free. This 
NetLogo implementation of the Kuramoto model sets up 
randomly selected initial phase angles and intrinsic fre- 
quencies drawn from a uniform distribution in the range 
0.02zLw/2. Thus the width w of the distribution, the cou- 
pling a, the number of nodes N and the type of graph 
(using p) are selected prior to execution of a simulation. 
The position of phases around the circle are indicated in 
the simulation graphic by a rotating colour scheme. In 
this model, "correlations" are plotted as a function of 
time based on the original order parameter of Kuramoto 
[4], which can easily be shown to display the same dy- 
namical behaviours of our r considered above. The model 
is designed so that for any graph that can be dialled up a 
full range of behaviours, from incoherence up to complete 
synchronisation, can be exhibited for a reasonable sweep 
of the parameters of coupling and frequency distribution 
width. 

We comment in more detail in the next section on the 
Laplacian spectral properties of various well-known com- 
plex networks. Here it suffices to say that for demonstrat- 
ing the transition from incoherence to synchronisation we 
have selected a scale-free graph {p = 5) with = 20 
nodes. The reason for this choice is to avoid high sup- 
pression by inverse powers of N of potentially oscillatory 
behaviour in r, see Eqs.( 27|28 l, while on the other hand 
avoiding a large build up of low-lying Laplacian modes. 
We thus have a system close to being consistent with the 
conditions of our approximations. 

A number of screen captures from simulations are pre- 
sented. We first show simulation for the lowest possi- 
ble coupling in the model settings, a = 10^^, in Figjl] 
Though r can occasionally reach large values there is no 
evidence of any order in the evolution of the fiuctuations 
which is consistent with incoherence. 

We next show two cases at coupling of cr = 2 x 10~* 
at which quite different periodic behaviours in the order 
paramater can be seen. In the first instance, Fig|2] there 
is evidently a single oscillating mode. This behaviour 
has been reproduced a number of times though each case 
shows slightly different profiles according to the (small) 
number of nodes participating in such modes. However, 
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FIG. 1: NetLogo simulation of the Kuramoto model for 
a scale-free network of 20 nodes. The coupling is set to 
a = 10~*. The colours of the nodes in the network dia- 
gram visually represent the phase angles of the oscillators 
and give a sense of how many oscillators are locked (the green 
nodes in this snapshot) versus how many are drifting (other 
colours). The order parameter "correlation" shows incoherent 
behaviour consistent with the poorly coupled network failing 
to synchronise. 



other instances, for example shown in Figj3] exhibit two 
superposed periodic modes. Again, a number of vari- 
ations on this are observed for different random initial 
conditions and intrinsic frequencies. 
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FIG. 2: NetLogo simulation of the Kuramoto model for 
a scale-free network of 20 nodes. The coupling is set to 
CT = 2 X 10~*. The order parameter "correlation" shows pe- 
riodic behaviour consistent with the limit cycle oscillation of 
a single mode. From the colours of the nodes we can see that 
many nodes are locked (green) while several (yellow and red) 
have not. This permits us to conclude that the limit cycle be- 
haviour is not a consequence of a single oscillating node but 
of collections (albeit in small numbers) of nodes. 

Finally, the coupling is increased to ct = 5 x 10^*^ 
whence rapid synchronisation occurs as seen in Fig|4j 

These behaviours are consistent with those identified 
analytically for the one and two mode cases so that 
the periodicity is rather robust though the precise form 
it takes can vary according to randomness within the 
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FIG. 3: NetLogo simulation of the Kuramoto model for a 
scale-free network of 20 nodes. The coupling is set to cr = 
2 X 10~*. The order parameter "correlation" shows periodic 
behaviour consistent with the limit cycle oscillation of two 
modes. 
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FIG. 4: NetLogo simulation of the Kuramoto model for a 
scale-free network of 20 nodes. The coupling is set to cr = 
5 X 10"*. The order parameter "correlation" shows rapid 
synchronisation. 



network and selection of intrinsic frequencies. Unfortu- 
nately, we are not able to extract the network properties 
directly from NetLogo in order to check the low-lying 
Laplacian spectral properties for these cases. 



VI. DISCUSSION 

We have considered second order approximations 
about the synchronised fixed point in the Kuramoto 
model on a general network. Motivated by insights from 
simulations that network sub-structures with smallest 
Laplacian eigenvalues approach synchronisation later in 
time, we have assumed higher modes to have reached 
steady state and analysed the dynamical equations for 
the remaining low- lying normal Laplacian modes. These 
equations resemble logistic and Lotka-Volterra equations 
and have yielded critical values for the coupling con- 
stant distinguishing behaviours such as convergence to 



the fixed point and limit cycles. In particular, this criti- 
cal coupling is inversely proportional to the square of the 
algebraic connectivity. For larger numbers of low-lying 
modes participating in the dynamics more exotic dynam- 
ics can be exhibited, such as tori and chaos. In support of 
these results for one and two modes we have used NetL- 
ogo simulations to show evidence of periodic behaviour as 
the coupling is increased. This shows that for dynamical 
processes on networks there can be intermediate regimes 
of behaviour, between order (synchronisation) and chaos 
(incoherence), exhibiting rich structured behaviour, an 
example of Edge of Chaos phenomena. 

At a broader level, these results show that for general 
networks there is no single critical coupling for the onset 
of synchronisation but many such coupling values indi- 
cating the regimes within which substructures can sta- 
bilise into states of partial synchronisation. We expect a 
similar hierarchy of couplings to be reflected in the non- 
perturbative regime, where the bulk of the normal modes 
have their dynamics "extinguished" . Our lack of analyti- 
cal handle on the mechanism at play in this regime means 
we cannot yet compare with the case where this model 
all began, the complete graph: the high degeneracy of 
Laplacian eigenvalues here (infinite degeneracy for Ku- 
ramoto's N ^ oo case) means that the normal modes' 
capacity to reach steady state lies entirely in this non- 
perturbative regime. 

However, we can speculate to some degree on the role 
of Laplacian modes here. We notice that, whereas for 
regular graphs the Laplacian spectra are quite dispersed 
with multiple gaps between groups of (often degenerate) 
eigenvalues, for all known complex networks there oc- 
curs a clear accumulation of modes about some eigen- 
value 28J: a "bulk" of modes is present in the spec- 
trum. Erdos-Renyi random graphs generated by a low 
probability of creating a random link between any two 
nodes show this bulk at small eigenvalues while for large 
probability this bulk shifts up to higher values. Watts- 
Strogatz networks with small probability for rewiring a 
lattice graph have the bulk in a narrow curve around the 
mean nodal degree with smaller peaks reflecting the orig- 
inal spectrum for the lattice while for larger probability 
the main peak broadens. Finally, scale-free graphs have 
the main bulk occuring narrowly around eigenvalues of 
order unity with a few eigenvalues less than one. This 
bulk in the Laplacian spectrum is the common feature of 
any graph which, according to numerical studies, read- 
ily demonstrates synchronisation at finite values of the 
coupling. 

A nonperturbative mechanism is therefore suggested: 
close-lying eigenmodes somehow mutually interfere with 
each other such as to cause them to de-excite to their 
steady states; when the majority of eigenvalues are 
closely spaced this de-excitation results in a large num- 
ber of modes reaching steady state values. The thresh- 
old for this de-excitation is given by a critical coupling 
value. For the complete graph, the complete degener- 
acy of modes means this de-excitation across modes is 
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nearly instantaneous. For other complex graphs this de- 
excitation takes some time to propagate, leaving some 
residual number of dynamical low-lying modes whose 
suppression has been the subject of this paper. 

This insight suggests new approximations of the dy- 
namical equations which may provide a handle on this 
nonperturbative mechanism. Our future work will seek 
to provide more substantial support to this hypothesis. 
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